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Abstract This paper presents a novel method for reti- 
nal vasculature extraction, using a vessel tracking method 
based on biologically inspired multi-orientation anal- 
ysis. We apply the multi-orientation analysis via so- 
called invertible orientation scores, modeling the corti- 
cal columns in the visual system of higher mammals. 
This allows us to generically deal with many hitherto 
complex problems inherent to vasculature tracking, such 
as tracking crossings, bifurcations, parallel vessels, ves- 
sels of varying widths and vessels with high curvature. 
The method runs fully automatically and provides a de- 
tailed model of the retinal vasculature, which is crucial 
as a sound basis for further quantitative analysis of the 
retina, especially in screening applications. 

Keywords Vessel tracking • Retinal imaging • 
Orientation scores • Retinal vasculature • Gabor 
wavelets 

1 Introduction 

1.1 Retinal vessel tracking 

The retinal vasculature is the only part of the body's 
circulatory system that can be observed non-invasively 
by optical means. A large variety of diseases affect the 
vasculature in a way that may cause geometrical and 
functional changes. Retinal images are therefore not 
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only suitable for investigation of ocular diseases such as 
glaucoma and age-related macular degeneration (AMD), 
but also for systemic diseases such as diabetes, hyper- 
tension and arteriosclerosis. This makes the retinal vas- 
culature a rewarding and well researched subject and 
and a growing number of image processing techniques 
are developed to segment and analyze the retinal vas- 
culature [TJ[2] 

Typically there are two types of methods for vessel 
extraction: pixel classification methods [3,4,5 and ves- 
sel tracking methods [6, 7,8,9, 10, lTJ[T2]. The first type 
of method classifies pixels as either being part of a ves- 
sel or background, resulting in a pixel map where white 
pixels represent blood vessels. Often a number of fea- 
tures are calculated for each pixel which are then used 
to classify, either supervised or unsupervised, the pixel 
as vessel or non-vessel. The other type of method, ves- 
sel tracking, is based on recursively expanding a model 
of the vasculature from a set of seed points. One advan- 
tage of vessel tracking over pixel classification is that it 
guarantees connectedness of vessel segments whereas in 
pixel classification methods this is not necessarily the 
case. For further quantitative analysis of the vascula- 
ture, tracking algorithms are preferred because they in- 
trinsically provide geometrical and topological informa- 
tion. For example, vessel widths, curvatures, segment 
lengths, bifurcation density and other features can rel- 
atively easily be extracted from the generated vessel 
models. Some of the most commonly investigated fea- 
tures of the vasculature are the arteriolar to venular 
ratio (AVR) [13 , curvature [14 (often incorrectly re- 
ferred to as tortuosity, tortuosity implies 3D geometry) 
and junction features such as bifurcation angle and the 
junctional exponent [15J[T6]. 

Quite a number of papers on vessel tracking have 
been proposed, each with different approaches. There 
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(a) Image 



(b) Orientation score 



(c) Image 




U f {x, y\ 0) 



(d) Orientation score 




(e) Transforms 



(f) Edge tracking in OS 



Fig. 1: (a,b): Crossing lines in the image domain (a) are clearly unwrapped in the orientation score domain (b). (c,d): The 
domain of an orientation score is curved and the orientation dimension is 2tt -periodic, (e) An invertible orientation score 
transformation neatly distributes pixel information from an image f over a set of orientations in an orientation score U f . 



that allows 



An orientation score transformation is invertible if there exists a well-posed inverse transformation 
exact reconstruction from the orientation score domain, (f) Disentanglement of crossing structures in the score domain allows 
vessel tracking through crossings. 



are methods based on active contours [6lfTU]. matched 
filters [7,9 , 12 j, probabilistic tracking [8] and many more 
[8 ,11. The majority of papers on vessel tracking report 
limitations regarding tracking blood vessels through cross- 
ings, bifurcations and/or more complex situations. As 
with so many visual tasks, the performance of comput- 
erized vessel tracking algorithms is still inferior to that 
of human observers. In order to deal with the complex 
situations that occur during vessel tracking, we pro- 
pose a novel vessel-tracking approach inspired by recent 
findings of the orientation processing in early stages of 
the visual system. More specifically, our vessel tracking 
method is based on invertible orientation scores, model- 
ing the cortical columns in the visual system of higher 
mammals. The method provides a detailed model of 
the retinal vasculature, where every branch is named 
and analyzed individually while keeping the hierarchi- 
cal structure. This model provides a good basis for fur- 
ther analyses of the retina. 



1.2 Brain-inspired orientation scores 

In 1959, Hubel and Wiesel [17] showed that the visual 
cortex in the mammalian brain contains so called sim- 



ple cells, whose receptive fields are tuned to various 
locations and orientations. They showed that the orien- 
tation preference of simple cells gradually varies across 
the surface of the primary visual cortex. Optical imag- 
ing techniques have shown that assemblies of oriented 
receptive fields are grouped together in pinwheel-like 
structures in the cortical columns [18], known as the 
orientation preference structure. The orientation pref- 
erence structure is a mapping of the Euclidean motion 
group space R 2 x T onto the 2D surface of the pri- 
mary visual cortex, where T denotes the orientation 
domain. Due to the difference in dimensionality, the 
orientation preference structure is punctuated by so- 
called pinwheels, which centers are singularities in this 
mapping [19 . 

Inspired by the orientation selectiveness of the brain, 
Duits et al. [20 developed a mathematical framework 
for image processing and analysis based on 2D orienta- 
tion scores^\ Similar to the perceptual organization of 
orientation in the visual cortex, a 2D orientation score 
is an object that maps 2D positions and orientation an- 

1 We use the name orientation score rather than an ori- 
entation space, as the latter term is inappropriate given the 
dimensionless orientation parameter. 
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gles (x, y, 6) to complex scalars. Instead of assigning an 
orientation to each position and thereby extending the 
codomain, we extend the domain of the image where 
our modeling can deal with multiple orientations per 
position. When constructing an orientation score it is 
crucial one does not tamper the data evidence before 
tracking takes place. Therefore we need invertible orien- 
tation scores that provide a full comprehensive overview 
of how the image is decomposed out of multiple orien- 
tations. In invertible orientation scores, all orientated 
structures are disentangled, see Fig. [I] Our new ves- 
sel tracking algorithm acts on the orientation score of 
an image, rather than on the image itself, an approach 
which overcomes many problems that occur during con- 
ventional vessel tracking at crossings. 

The article is structured as follows: First, in sec- 
tion [2j theory about orientation scores is provided. In 
section [3J two vessel tracking approaches based on ori- 
entation scores are described: 

— the ETOS-algorithm: an all-scale approach based on 
a new class of wavelets, the so-called cake wavelets 

— the CTOS-algorithm: a multi-scale approach based 
on the classical Gabor wavelets 

The ETOS algorithm performs best with invertible ori- 
entation scores based on all-scale cake wavelets [21, 
[22] (in comparison to non- invertible orientation scores 
based on Gabor wavelets). The second approach re- 
quires a multi-scale and orientation decomposition. Both 
approaches are described in section pTH and subsequently 
in section |3.2| they are evaluated. It will turn out that 
ETOS based on cake wavelets has several advantages 
over CTOS based on Gabor wavelets. We have vali- 
dated ETOS more extensively by comparing it to the 
state of the art in retinal vessel tracking [6 ,23 ,24 using 
the publicly availabe REVIEW database [6] . For precise 
quantitative analysis we must capture the complete vas- 
culature, a process which we call vasculature tracking. 



In section \4~T\ we describe our vasculature tracking al- 
gorithm, composed of proper initialization, junction de- 
tection and junction resolver algorithms. In section |3~2 
the correctness of topology of the models is evaluated 
using images of the HRFI-database [25]. General con- 
clusions can be found in section [5] 



2 Orientation scores 

Orientation detection and encoding is a common sub- 
ject in image processing. E.g., Frangi et al. [26] detect 
retinal blood vessels by calculating a vesselness value, 
obtained by eigenvalue analysis of the Hessian matrix. 
In [27] this is done by taking the maximum modulus 



over a set of oriented Gabor wavelet responses. Be- 
sides the presence of local orientations, the orientation 
value (s) may be relevant as well. For instance in ves- 
sel tracking and certain other image enhancement tech- 
niques like coherence enhanced diffusion [28] and orien- 
tation channel smoothing [29] 

The most commonly used methods to detect orien- 
tations are capable of detecting only one orientation per 
position. However, in practice a position may contain 
multiple oriented structures, e.g. at corner points, cross- 
ings and bifurcations. By using oriented wavelets and 
steerable filters [30 ,31,32 ,33 orientation confidence mea- 
sures can be extracted for any given orientation, thus 
allowing for the detection of multiple orientations per 
position. Oriented wavelets allow for a transformation 
from an image to an orientation score, where each com- 
bination of position and orientation is mapped to a sin- 
gle value [3l[20][22], see Fig. [I] 

In his pioneering paper, Kalitzin [34] proposed a 
specific wavelet, given by 



t/>(x,y) = 



N 

£ 

n=l 



2 * = 



x + iy, 



that, by approximation, guaranteed invert ibility from 
the orientation scores back to the original image, with- 
out loss of information. These wavelets belong to a spe- 
cific class of so called proper wavelets (wavelets that 
allow well-posed reconstruction [20]|22]), and are found 
by expansion in eigenfunctions of the harmonic oscilla- 
tor. The advantage of this expansion is that this steer- 
able basis is Fourier Invariant, allowing to control the 
wavelet shape in both the spatial and Fourier domain. 
The disadvantages of such kernels are however that 1) 
their series do not converge in L2 and truncation of 
the pointwise converging series heavily affects the shape 
and induces undesirable oscillations [2TJ p. 140-142]; 2) 
the wavelet explodes with ~ (87r) 1 / 4 v / r(l — 0(r~ 2 )) in 
the radial direction along its orientation [2TJ App.7.3]; 
3) does not allow approximate reconstruction by inte- 
gration over S 1 . 

In Kalitzin's paper the well-posedness of the re- 
construction was not quantified. Analysis of the well- 
posedness was done by Duits using the function M^, 
which will be explained in section 2.1. The function 
M-0 indicates how well spatial frequencies in an im- 
age are preserved after a transformation, and it can be 
seen as a measure for stability of the inverse transfor- 
mation [22 ,35,21,20 . Within these papers a new class 
of proper wavelets called cake wavelets are presented. 
These are oriented wavelets, able to capture all image 
scales without any bias to a specific scale. This property 
is crucial in our vessel edge tracking approach since it 
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eliminates the need for a multi-scale approach, as will 
become more clear in the next sections. 

In contrast, the Gabor wavelet is another oriented 
wavelet, which is widely used in the field of image pro- 
cessing because of its capability to detect oriented fea- 
tures at a certain scale. An example of segmenting cross- 
ing lines using Gabor wavelet based orientation scores 
is given in [36] . The property to tune the Gabor wavelet 
to capture features at a specific scale can be very use- 
ful, but the scale selection also implies exclusion of other 
scales. The single-scale Gabor transformation causes in- 
formation from the original image to be lost and the 
transformation is therefore non-invertible. In order to 
introduce invertibility, one has to use a multi-scale ap- 
proach [37] , which is also computationally more expen- 
sive. 

Since no information should be lost during the orien- 
tation score transformation : / — » Uf (see Fig. le), 
the notion of invertibility of an orientation score trans- 
formation is essential. The invertibility allows us to re- 
late operators on images to operators on orientation 
scores, and vice versa. Using invertible orientation scores, 
one can employ the automatic disentanglement of lo- 
cal orientations involved in a crossing (cf. Fig. [Tf|). For 
line/contour enhancement, this has lead to a generic 
crossing preserving diffusion method [38,39,40 which 
outperformed related diffusions acting directly on the 
image domain. For tracking of crossing blood vessels a 
similar advantage can be employed, as we will show in 
this article. 



2.1 Construction of orientation scores 



Consider a 2D image / as a function / : R 2 — >• R, with 
compact support on the image domain Q = [0,X] x 
[0,y], with image dimensions X, Y G R + , and which 
we assume to be square integrable, i.e. / G L2(R 2 ). An 
orientation score, constructed from image /, is defined 
as a function Uf : I 2 x S 1 — >> C and depends on two 
variables (x, 0), where x = (xi,^) £ ^ 2 denotes posi- 
tion and 6 G [0, 2tt] denotes the orientation variable. 

An orientation score Uf := W^f of a function / 
can be constructed by means of convolution with some 
anisotropic wavelet ip via 



m n ( P 2 r 1 ) 



l/ / (x,fl) = (^*/)(x) 



JR 2 



R e - 1 (y-x))/(y)rfy, 



(1) 



where ip G L2(R 2 ) is the convolution kernel with orien- 
tation = 0, i.e. aligned with the vertical axis in our 
convention, and where denotes the transformation 




Fig. 2: Plots of Mn {p 2 /t), with t = - \' ^ r for N = 5, 10, 



15, 20, 25 



1 + 2AT 



between image / and orientation score Uf. The over- 
line denotes complex conjugate, ^(x) = ^e{~ x ) an d 
the rotation matrix R<9 is given by 



R/Q 



cos — sin 



sine 



cost 



(2) 



see Fig. le Exact reconstruction^ from the orientation 
scores constructed by ([!]) is given by 



(3) 



where T is the unitary Fourier transform on R 2 , where 
denotes the adjoint wavelet transformation (see [21 



for details), and M ( 



— » R + is calculated by 



2tt 



2ir 



me]\ 2 d0. (4) 



The function provides a stability measure of the 
inverse transformation. Theoretically, reconstruction is 
well-posed, as long as < 5 < M^(uj) < M < oo where 
6 is arbitrarily small, since then the condition number 
of is bounded by M 5~ x [21j Thm. 20]. In prac- 
tice, to prevent numerical problems, it is best to aim 
at M^j(uj) ^ 1 for < £>, where g is the Nyquist fre- 
quency of the discretely sampled image, meaning that 
all relevant frequency components within a ball of ra- 
dius g are preserved in the same way. Because of the 
discontinuity at ||u;|| = g, which causes practical prob- 
lems with the discrete inverse Fourier transform, we 

2 The reconstruction formula can easily be verified using 
the convolution theorem, J 7 [f*g] = — •^ r [/]^ r [p] 5 an d the fact 
that T 
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will focus on wavelets -0, with M^{ui) = Mn (p 2 t x ) 
N G N, t > and p = \\lj\\ where 



M N ( P 2 t- i ) = e t j2 [p ' < t 



fe=0 



(5) 



where £ denotes a scale parameter. To fix the inflec 
tion point close to the Nyquist frequency, say at p - 

2( 7 p) 2 



7£> with <C 7 < 1, we set t 



1 + 27V 



-j-^M N (p 2 t 1 )\ p=1Q = 0). The function Mn basically 

is a Gaussian function at scale t, multiplied with the 
Taylor series of its inverse up to a finite order 2N to 
ensure a slower decay. The function M n smoothly ap- 
proximates 1 on the domain p G [0, g], see Figure [2] A 
wavelet ip : R 2 —> C with such a will be called a 
proper wavelet. 

The methods presented in this paper are best de- 
scribed by a moving frame of reference that lives in the 
domain R 2 x S 1 of an orientation score: 



= — sm(6)e x + cos(0)e y 
= cos^e^ + sm(0)e y 
e* = (0,0,1), 



(-sin(0), cos((9),0) 
(cos(<9),sin(<9),0) 



(6) 



with r] = — s'm(6)x + co$(9)y, £ = cos(6)x + sm(9)y, 
and with e x = (1, 0, 0) and e y = (0, 1, 0), such that at a 
given point in the orientation score (x, 2/, 0), points in 
the direction of and perpendicular to it, see Fig.[Tc| 
and Id We will often rely on the notation in Eq. ^ in 



a rigid body motion, SE(2) is equipped with the group 
product: 



g-sf = (x,R^)(x / ,R^) = (lW + x,IW), 



(8) 



which is consistent with Eq. The moving frame of 
reference (|6| corresponds to the so-called left-invariant 
vector fileds in SE(2). For details see [40j Fig. 2.5a]. 



(so that 2.3 Cake wavelets 



Cake wavelets are constructed from the Fourier domain. 
By using polar coordinates, the Fourier domain can be 
uniformly divided into N Q equidistant and equally wide 
samples ("pieces of cake") in the angular direction. The 
spatial wavelet is given by 



^ cofce (x)=^- 1 [^ cafce ](x)G aa (x), 



(9) 



the remainder of this article. 



where G as is a spatial Gaussian window, with < 1 << 
cr s , that is used to avoid long tails in the spatial do- 
main. Note that multiplication with a large window in 
the spatial domain corresponds to a convolution with 
a small window in the Fourier domain, such that 
is hardly affected with a s sufficiently large. Function 
ip cake is given by 

~ , ( (u> mod 2tt) - tt/2\ 

with uj = (p cos (/?, p sin ip) and where sq = 2ttN~ 1 is the 
angular resolution in radians. The function M n speci- 
fies the radial function in the Fourier domain given by 
([5]). Bk denotes the kih. order B-spline given by 



2.2 The domain of orientation scores: SE(2) 

The domain of an orientation score is the set R 2 x S 1 . 
However, from Fig. [id] one can recognize a curved ge- 
ometry on the domain of orientation scores. This is re- 
flected in the fact that and vary with 0. Mathe- 
matically, this is modeled by imposing a group structure 
on the set R 2 x S 1 . This group structure comes from 
rigid body motions g = (x, R^) acting on R 2 x S 1 via 

0. (x',0') = (R*x' + x,0 + 0')- (7) 

Note that (x, 0) = (x,R<9)(0,0) which allows us to 
uniquely identify 

R 2 x S 1 3 (x, 0) o (x, He) = 9 e R 2 x 50(2), 

1. e. to identify the space of positions and orientations 
with the rigid body motion group SE(2) = R 2 x SO(2). 
As the combination of two rigid body motions is again 



B k (x) = (B k _ 1 *B )(x), 

f 1 if -1/2 < * < +1/2 ^ (11) 
\ otherwise 

Orientation scores constructed from an image / using 
cake wavelets are denoted by U c ^ ake . 

The approach of constructing wavelets directly from 
the Fourier domain allows indirect control over the spa- 
tial shape of the filter, and it can easily be adapted. For 
example, the number of orientations N Q specifies the an- 
gular resolution sq: If N Q is large, the resolution in the 
orientation dimension is large and the filters become 
very narrow. This is illustrated in Figure |4j Moreover, 
because B-splines and the function M n are used to sam- 
ple the Fourier domain, the sum of all cake wavelets is 
approximately 1 over the entire Fourier domain (Fig. [3] 
and|4| and thus they indeed constitute proper wavelets, 
allowing stable reconstruction via Eq. In Eq. ^ 
one can omit division by M^ 1 « 1 in which case sta- 
ble reconstruction is obtained both by integration over 



6 



Erik Bekkers et al. 




Fig. 4: Overview of the wavelets used in this paper. From left to right: the real and imaginary parts of the wavelet in the spatial 
domain (zoomed by a factor of 8 for the scae of visualization), the wavelet in the Fourier domain and the an illustration of 
the Fourier domain coverage by the filters where contours are drawn at 80% of the filter maximum. Note that this last figure 
also gives an impression of M^. The top row shows the cake wavelet constructed using N Q = 36, middle row with N Q = 12 
and the bottom row shows the Gabor wavelet at scale a = 6/tt and N Q — 36. 



SE(2) := R 2 y\ S 1 and/or its discrete subgroup '. 



n(2iri/N) | 



0,1,..., TV - 1}. The cut-off 



with Slf 

frequency (the inflection point) of the function M ni 
which is usually set as the Nyquist-frequency, could be 
lowered to filter out high-frequency noise components. 

Cake wavelets are quadrature filters, meaning that 
the real part contains information about the locally 
even (symmetric) structures, e.g. ridges, and the imag- 
inary part contains information about the locally odd 
(antisymmetric) structures, e.g. edges. That is, the real 
and imaginary part of the filter iJjq are related to each 
other by the Hilbert transform in the direction perpen- 
dicular to the wavelets orientation, which is defined by 

H"(^)(x) = ^- 1 [a;^i sign(^ e?) )J^](u;)](x), (12) 

where specifies the direction in which the Hilbert 
transform is performed, recall Eq. (TOT) . 



The quadrature property is useful in our vessel track- 
ing approach, since it allows us to directly detect vessel 
edges from the imaginary part of the orientation scores, 
without having to calculate first-order derivatives per- 
pendicular to the vessel orientation. In our applications 
we did remove the DC-component for the real part to 
avoid responses on locally constant images. Figure [4] 
shows the real and imaginary part of the cake wavelet 
in the spatial domain, as well as the coverage of the 
wavelet in the Fourier domain. 

A final remark on cake wavelets: Since the cake 
wavelets uniformly cover the Fourier domain {Mn ~ 
1), they allow us to use a fast approximate reconstruc- 
tion scheme, which is given by integration of the orien- 
tation scores over the angles only: 



r 



c (x) 



i 

2^ 



m ake (x,e)d0. 



(13) 
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Fig. 3: TTie lise o/ B-splines in the construction of cake 
wavelets. Plot showing quadratic B-splines (k=2), the sum 
of all shifted B-splines add up to 1. The image in the up- 
per right corner illustrates a Fourier cake wavelet ip cake (uj) 
constructed using qua dratic B-splines and Mn with N = 60, 
according to Eq. 



2.4 Gabor wavelets 

In the previous section we discussed the cake wavelet, 
an all-scale oriented wavelet which can be used to con- 
struct invertible orientation scores. In this section we 
will discuss the Gabor wavelet, a scale selective oriented 
wavelet, which can be used to construct non-invertible 
orientation scores. The Gabor wavelets are directional 
wavelets, which can be tuned to specific spatial frequen- 
cies (and inherently scales). In the field of retinal image 
processing they are used for vessel detection in various 
studies [27,41 . We can exploit the tuning of the wavelet 
to specific spatial frequencies to match differently sized 
blood vessels. The 2D Gabor wavelet is a Gaussian, 
modulated by a complex exponential (a wave function), 
and is defined as: 



^ Ga6or (x) 



1 

C-0 



e ik x e 2 



|Ax| 2 



(14) 



C^jj = 27Ty/e e 2 



lA^kol 2 



where A = diag\e~ x ' 2 ^ 1] with e > 1 is a 2 x 2 diagonal 
matrix that defines the anisotropy of the wavelet. The 
vector ko defines the spatial frequency of the complex 
exponential and normalizes the wavelet to unity. 
In our method we use e = 4, which makes the filter 
elongated in the x-direction and we choose kg = (0,3), 
which causes oscillations perpendicular to the orienta- 
tion of the wavelet. 

We can dilate the filter by a scaling parameter a > 0: 



Orientation scores constructed from an image / using 
Gabor wavelets at scale a are denoted by Uf£ bor . 

The real and imaginary part of the wavelet in the 
spatial domain, as well as the coverage of the wavelet in 
the Fourier domain are shown in Figure [4j Similar to the 
cake wavelets, Gabor wavelets also have the quadrature 
property orthogonal to their orientation. 

In the Fourier domain, the Gabor filters are rep- 
resented as Gaussian functions shifted from the origin 
by ko. The set of all rotated Gabor functions at a cer- 
tain scale a covers therefore a certain annulus in the 
Fourier domain, and a single Gabor wavelet can thus 
be regarded as an oriented band-pass filter. This is also 
clearly depicted by the outlines of the frequency re- 
sponses as shown in Figure [4j 



2.5 Double-sided vs single-sided wavelets, orientation 
vs direction 

The cake and Gabor wavelets are double-sided wavelets 
which do not distinguish between a forward or back- 
ward direction (they are symmetric with respect to the 
y-axis). In order to distinguish between tt symmetries 
and 2tt symmetries (see Fig. [5|, and to be able to han- 
dle bifurcations, we decompose the orientation scores 
into a forward and backward direction^] denoted by a 
+ and — symbol respectively: 



U f (x, y, 0) = U-f{x, y, 0) + UJ (x, y, 0), 
where 



C/+(x, 9) = / R2 ^+(V(y-x)) /(y)dy, 
Uj{n, 9) = J R2 V(y " *))/(y)<*y, 

and where 

^+(x,y) = w(x)ip(x,y), 

i/)~(x,y) = w(-x)r/>(x,y) = (1 - w{x))ijj{x,y) 



with 



w(x) 



1 



-erf(x) 



l r x 

2Wo 



~y dy. 



(16) 



(17) 



(18) 



(19) 



2 v 7 2 2tt 

Note that by using the error function, we have ip = 
ip~ + ip + and UJ {x,y,6) = U^~(x, y, + 7r), so that 

Uf(x,y,0) = U~^(x,y,6) + Uf(x,y,0 + ir). It is thus 
possible to choose one of the single-sided wavelets to 
construct a directional orientation score, while still be- 
ing able to access the original (double-sided) orientation 



3 I.e. we extend the domain SE{2) 
orientation scores to the group E(2) 



a6or (x) = aT V Gabor (a _1 x). 



(15) 



^ 2 x SO(2) of our 
^ 2 x 0(2), where 
0(2) = {M e R 2x2 \M T = M' 1 } also includes, besides ro- 
tations (with det M = +1), reflections (with det M = —1). 
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(b) Single-sided 

Fig. 5: Comparison between double- and single-sided cake 
wavelets by visualisation of the orientation column £7/(x, •) 
at several points in a fundus image. The orientation column 
at a certain point x is visualized by drawing 36 lines, each at a 
certain angle < 6 < 2tt, and of which the length in direction 
is biven by the absolute value of the score \Uf(x.,6)\. 



score. Figure [5] demonstrates the advantage of using 
single-sided wavelets over double-sided wavelets in the 
case of direction estimation based on the orientation 
column of a score Uf(x,y,-). 



Ae, 



le f ,k-l 



(V/, 



Ck-1 1 



e ?7,k-l 



\Ck-l 



Fig. 6: Schematic overview image of ETOS. Using the de- 
tected vessel center point Cfc_i and the orientation 6k-i de- 
tected at the vessel edges Ufc_i and v^-i at iteration k — 1, 
a rough estimation of the next center point found by step- 
ping forward in the vessel direction e^^-i with step size X. 
Through the estimated center point a line is defined on 
which the new vessel edges u/c and v/c are detected. At these 
edges the vessel orientation Ok is detected and the final pre- 
cise vessel center point is calculated as the mean of the 
two edges. 



3 Vessel tracking using orientation scores 

This section describes two methods for vessel tracking 
based on orientation scores constructed by the wavelets 
described in section [2j 

1. The ETOS method: Edge Tracking based on Ori- 
entation Scores. The method is based on tracking 
vessel edges through the orientation scores and is de- 
scribed in section 13.1.11 In section [3.2. II the method 
is qualitatively evaluated; Here the method's per- 
formance is tested on orientation scores generated 
by cake wavelets as well as by Gabor wavelets to 
compare the difference between invertible and non- 
invertible orientation scores respectively. In section 
|3.2.2[ the edge detection capability of the method is 
quantitatively compared to state-of-the-art retinal 
vessel width measurement algorithms. 

2. The CTOS method: Centerline Tracking based on 
multi-scale Orientation Scores. This method is spe- 
cially designed to work with multi-scale orientation 
scores and exploits the scale selective properties of 
the Gabor wavelets to filter out a vessel's central 



light reflex. The method is described in section 3.1.2 
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and qualitatively evaluated and compared with ETOS 
in section 13.2. II 



3.1 Methods 

3.1.1 ETOS: Edge tracking using orientation scores 

In orientation scores information from the original im- 
age is neatly organized in different orientations, lead- 
ing to the disentanglement ("pulling apart") of crossing 
structures. Moreover, because of the quadrature prop- 
erty of the wavelets used in the construction of orien- 
tation scores, important edge information is well repre- 
sented in the imaginary part of the score. We therefore 
propose the ETOS algorithm, which directly acts on 
the orientation score, and which is based on tracking 
the vessel edges. 

The ETOS algorithm iteratively expands a blood 
vessel model by detecting, at each forward step fc, the 
optimal edge locations g Uk = (u k ,6 k ),g Vk = (v k ,9 k ) e 
SE(2) from the orientation score, where and v/e de- 
note the 2D left and right edge position respectively, k 
denotes the orientation of the blood vessel, and g Uk and 
g Vk are the group elements representing respectively the 
left and right edge in the orientation score. At each it- 
eration the vessel center point c k and the vessel width 
w k are defined as follows: 



W k = 1 1 U]fe - Vfc||. 



(20) 
(21) 



To describe our method we will rely on a moving frame 
of reference with basis vectors e^ fc , e^ fc and e^ fe , which 
are described by the orientation parameter k and Eq. (J6|. 

In our method the edge positions g Uk and g Vk are de- 
tected in the orientation score from the tangent plane 
(plane tangent to the vessel orientation) spanned by 



eg and (yellow plane in Fig. 8a). An edge can be 
detected as a local optimum from the imaginary part 



of this plane (see Fig. 8c); a local minimum and max- 



imum for the left and right edge respectively. Such a 
local optimum basically describes where the wavelet 
used to construct the orientation score is best aligned 
with the vessel edge (in terms of position and orienta- 
tion). More information on optimization processes and 
optimal paths in orientation scores can be found in ap- 
pendix [A| 

For the sake of speed and simplicity, the process of 
detecting the optimal edge positions is separated into 
two ID optimization tasks which simply involve the de- 
tection of local minima (left edges) and maxima (right 
edge); first the edge locations and are optimized 



in the direction, then the corresponding orientation 
is optimized in the direction (Fig. [8}>d). By con- 
sidering the continuous properties of the blood vessels 
(e.g. continuous vessel widths), we use a paired edge 
tracking approach where the left and right edges are 
detected simultaneously. This approach has the advan- 
tage that even if one of the edges is less apparent in the 
image (e.g. at crossing points, parallel vessels and bifur- 
cations), both edges and their orientation can still be 
tracked. A weak point of this approach is however that 
abrupt changes in vessel width (e.g. at stenoses and 
aneurysms) may become unnoticed or detected with 
less detail. 

Based on a-priori knowledge about the previous ves- 
sel orientation, edges and center point, as stored in 
0fc_i, Ufc_i, Vfe_i and c k -i respectively, phase 1 is started 
by estimating a new vessel center point c k according to 



c k = Cfc-i + A e^ fc _ 1 



(22) 



where A (typically in the order of 2 pixels) is the track- 
ing step-size. New edge points are selected from a set 
of points Pk(j]) on a l me L k going through the esti- 
mated center point and perpendicular to the vessel 
orientation 6 k -\\ 



L k = {Pk(v) I V ^ [-VmaxiVmax]}, 

with 

Pk(v) = Cfe +r/e r?fc _ 1 , 



(23) 



(24) 



where 77 is a parameter describing the distance to the 
estimated vessel center point and r] max > \\u k — v k \\ 
denotes the maximum distance to the estimated the 
vessel center point. Note that orientation k -i of the 
previous iteration is used as the new orientation is yet 
to be detected. 

An intensity profile I k (rj) can then be obtained from 
the orientation scores according to 



rHiv) = C//(p fe (r?) A-i). 



(25) 



An example of such an intensity profile is shown in Fig- 
ure 8d New edge points can now easily be found by de- 



tecting the local optima on the imaginary part of this 
profile. However, the detection of vessel edges is made 
more robust by taking into account that the vessel wall 
is a continuous structure, and that the width of a blood 
vessel gradually changes, rather than abruptly. There- 
fore, for the robust detection of edges from the inten- 
sity profile, we introduce the edge probability envelope. 
The edge probability envelope is used to indicate the 
most likely position of the vessel edges and it consists 
of two Gaussian distributions, one around the expected 
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(a) Image selection (b) Zoomed image /(■) (c) Uf ake (-,6 V 



(d) uf; a b °r(-,t 



(e) uf^{-,e v ) 



Fig. 7: Parallel blood vessels and orientation scores, (a) A selection of a fundus image and (b) a close-up view, (c-d) Slices of 
orientation scores constructed from (b) using cake wavelets and Gabor wavelets at scale a\ = 3 * 10/(27r) and a 3 — 3*30/(27r) 
respectively. The slices correspond to the orientation of the two parallel blood vessels, V . 




Orientation detection, I®(9) 
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(a) Vessels in OS and detection plane 



(b) Real 



(c) Imaginary 



(d) Detection profiles 



Fig. 8: Edge tracking in a ^-periodic orientation score constructed from double-sided wavelets, (a) Graphical representation of 
blood vessels in the orientation score. The real and imaginary part of the orientation score on the yellow plane V (perpendicular 
to the blood vessel) are represented in (b) and (c) respectively. In (c) the left and right edge of the blood vessel are expressed 
as black and white blobs respectively. The edge and orientation detection profiles are demonstrated in (d). 



left vessel edge position and one around the right vessel 
edge position. The envelope function is given by: 



w k w k 

Ew k ,*,>no(v) = -G a (v+ -Y-Vo) + G a (r]- — - r] ) 

lx 2 



G a (x) = 



1 



aV27r 



(26) 



where rjo is the estimated location of the vessel center 
on L/c, w k is the mean vessel width calculated over the 
last 10 iterations: 



10 



10 



1 10 

To ^ 

rn=l 



■v fc _ m ||, (27) 



the standard deviation of the Gaussian distributions 
is denoted with a, and which is used to align the 
envelope with the actual vessel profile, is an accurate 
estimation of the vessel center position on the scan-line. 
A typical edge probability envelope is given in Fig. [9b. 



An accurate value of 770 can be found by correlating the 
edge probability envelope E Wj(7jrio= o(r]) with the imag- 
inary part of the intensity profile and detecting 
a local maximum of this result. This is illustrated in 
Fig. [9J3. The left and right edges are finally detected 
as the arguments corresponding to the minimum and 
maximum points, Vieft anc ^ V right respectively, of the 
product of the envelope and the intensity profile: 

Vieft = argmin {Im I^rj) l&w^^ (r])\ } , 

rjE[-r)max,rjo] /og) 

Vright = argmax {Im ^(77) \E Wk ^ Vo (r])\ } , 

r}e[r l o,r lmax } 

see Fig. [9]i for an example. The new edge points can 
then be assigned by 



: Pk(Vleft)i 
PkiVright)' 



(29) 



Phase 2 now starts. An orientation Ok can be estimated 
by selecting the orientation that provides the highest 
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Fig. 9: Edge detection using the edge probability envelope, (a) Cross- sectional intensity profile taken from Fig. \ l0c\ showing 
many potential candidate left (L) and right (R) edge positions, (b) The edge probability profile, (c) Centering of the edge 
probability profile on the vessel of interest by means of correlation, (d) Enveloping the intensity profile results in clearly 
detectable left and right edge points. 
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Fig. 10: Scale selective orientation scores can filter out a vessels central light reflex, (a) A small sub-image showing a vessel 
with central light reflex and (b) the corresponding intensity profile, (c) A slice, corresponding to the vessel orientation V , of 
the orientation score constructed from (a) and (d) the corresponding intensity profile taken hereof. Note that from (d) the 
vessels center point can be roughly detected as a local minimum. 



orientation score response at the vessel edges. The ori- 
entation score response is a combination of the orienta- 
tion columns at the left and right edges. An orientation 
column at a certain position x is denned as E//(x, •). 
Since the imaginary part of the orientation score at a 
left edge is negative, and at a right edge is positive, the 
new vessel orientation is calculated as follows: 



% = argmax Im( 

6>G[0,2tt] 



-Uf{u k ,0) + U f (v k ,0)) 



(30) 



Finally the new center point c^, which may not be equal 
to Cfc, is calculated as the point between the two edges, 



according to Eq. (20). 



3.1.2 CTOS: Multi-scale vessel center-line tracking 
using orientation scores 



As described in section |2.4[ with Gabor wavelets it is 
possible to construct scale-specific orientation scores. 
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In this section the scale-selective property of the Gabor 
wavelets is exploited in the design of a fast orientation 
score based method called CTOS. The potential pres- 
ence of a central light reflex in a blood vessel makes 
the design of a simple and fast centerline tracking al- 
gorithm based on local minima tracking in the image 
nearly impossible. However, by using Gabor wavelets 
and appropriate scale selection, central light reflexes 
can be filtered out such that vessel center points can 
easily be detected as local minima on detection profiles 



(see Fig. [10]). 

The CTOS algorithm has a similar approach as the 
ETOS algorithm, but now each iteration k consists of 
three phases: A center point detection phase, an ori- 
entation Ok detection phase and an additional scale 
detection phase. Using the vessel center point Cfc_i, ori- 
entation Ok-i and scale a^-i, which were detected dur- 
ing iteration k — 1, phase 1 at iteration k is started by 



estimating the new center point as given by Eq. (22 ). 



The new center point is selected from a set of can- 



didate points Pk(v) as given by Eq. (24). From the 
candidate points Pk(v) a center point detection profile 



jV, Gabor ^ ^ g Qj^amed \yy. 



T ri,Gabor / \ 



(31) 



with Uf£^ i the orientation score generated by the Ga- 
bor wavelets at scale ctk-i- The new center point is 
detected as the coordinate belonging to the local mini- 
mum on the intensity profile nearest to c^: 



Cfe = Pfcfao), 
r]o = argmin 

max I 



Re I, 



r],Gabor 



(v)} 



In phase 2 the new vessel orientation is detected the 
local maximum of the negative orientation response, 




Fig. 11: A small section of a fundus image around the optic 
disk. Left: un-processed image. Right: high-pass filtered im- 
age obtained by estimating the local DC-value by Gaussian 
blurring the image with scale o — 32 and extracting it from 
the original image. 



nearest to the previous vessel orientation Ok-i' 

9 k = argmaxRe( -Uff k °\ (c k ,6) ). 
0e[o,27r] 



(33) 



Finally, at the new center point c& and orientation 
Ok , scale dk is detected as the scale that gives the largest 
negative response: 



a k = argmaxRe( -Uff or (c fc , k ) ). 



(34) 



Compared to the ETOS algorithm, this algorithm is 
very fast since it only requires a few simple (determinis- 
tic) detecting phases. However the combination of scale 
and orientation detection makes the algorithm slightly 
less stable: orientation detection depends on the correct 
detection of scale and vice versa. 



3.2 Validation 

The algorithms were tested on color fundus images. 
Processing of each image is done on the green channel, 
as this channel provides best contrast between blood 
vessels and background. For each image the luminos- 
ity is normalized by disposing low frequency luminos- 
ity drifts, a process which we will refer to as high- 



pass filtering (see Fig. 11). The low frequency drifts 
are detected by Gaussian blurring the image (typically 
a = 32), and are subtracted from the original image. 

3.2.1 Algorithm behavior at complex vessel junction 
points 



(32) A qualitative validation is done using a challenging set 



of 4 sub-images (see Fig. 12), which were taken from 



the high-resolution fundus images of the HRFI database 
[25] . This set of sub-images contains crossings, overlap- 
ping bifurcations with crossings, small vessels crossing 
large vessels, small vessels, curved vessels, parallel ves- 
sels, etc. The extracted set is therefore very suited to 
test the ability of the proposed tracking methods to deal 
with complex situations. In each sub-image we manu- 
ally placed seed points at the start of each blood vessel 
and at each bifurcation. In total 27 seed points were 
marked. Each seed point contains initial vessel center 
position, left edge position, right edge position and ori- 
entation, denoted by Co, uo, vo and 0q. The initial scale 
for the CTOS algorithm is detected as the scale that 
provides the largest scale response at cq and 0q ( see 
Eq. ([34])). 

The tracking experiments are conducted using the 
following set of tracking parameters: The step size is 
set to A = 2 pixels; The width of the scan line is set 
to 2r] rnax = 40 pixels; The number of orientations used 




Fig. 12: Results of vessel tracking on the test image set. From left to right: seed points, tracking results using the ETOS 
algorithm using invertible orientation scores (cake wavelets), tracking results of the ETOS algorithm using non-invertible 
orientation scores ( Gabor wavelets at scale r — 10 j and tracking results of the CTOS algorithm using orientation scores 
constructed at scales r — 5, 10, 15, 20, 25 and 30. Note that the results of the CTOS algorithm are only represented as 
centerlines since vessel width is not measured. From top to bottom: results on test image 1, 2, 3 and 4- 



to construct the orientation scores is set to N Q = 36 
and the standard deviation of the Gaussian distribution 
used for calculating the edge probability envelope is set 
to a = 3. 

The ETOS algorithm was tested on both invert- 
ible orientation scores, which were constructed by cake 



wavelets, and non-invertible orientation scores, which 
were constructed by Gabor wavelets. The scale of the 
Gabor wavelets was chosen in such a way that the rele- 
vant vessel features were presented as good as possible 
in the orientation scores (e.g. a scale too large would 
only represent the very large blood vessels correctly, 
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Table 1: REVIEW database comparison of successful measurement percentages (%), mean vessel widths (Mean) and standard 
deviations of the measurement errors (cr x ). 
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Standard 


100.0 


7.51 


0.00 


100.0 


13.79 


0.00 


100.0 


8.83 


0.00 


100.0 


4.35 


0.00 


Ol 


100.0 


7.00 


0.23 


100.0 


13.19 


0.57 


100.0 


8.50 


0.54 


100.0 


4.12 


0.27 


02 


100.0 


7.60 


0.21 


100.0 


13.69 


0.70 


100.0 


8.91 


0.62 


100.0 


4.35 


0.28 


03 


100.0 


7.97 


0.23 


100.0 


14.52 


0.57 


100.0 


9.15 


0.67 


100.0 


4.58 


0.30 


Gregson 


100.0 


7.29 


0.60 


100.0 


12.80 


2.84 


100.0 


10.07 


1.49 


100.0 


7.64 


1.48 


HHFW 


96.3 


6.47 


0.39 


0.0 






78.4 


7.94 


0.88 


88.3 


4.97 


0.93 


IDG 


100.0 


4.95 


0.40 


98.6 


6.30 


4.14 


99.9 


5.78 


2.11 


99.6 


3.81 


0.90 


2DG 


100.0 


5.87 


0.34 


26.7 


7.00 


6.02 


77.2 


6.59 


1.33 


98.9 


4.18 


0.70 


ESP 


100.0 


6.56 


0.33 


93.0 


15.70 


1.47 


99.6 


8.80 


0.77 


99.7 


4.63 


0.42 


Graph 


99.4 


6.38 


0.67 


94.1 


14.05 


1.78 


96.0 


8.35 


1.43 


100.0 


4.56 


0.57 


ARIA 


100.0 


6.30 


0.29 


100.0 


14.27 


0.95 


99.0 


8.07 


0.95 


99.5 


4.66 


0.32 


ETOS 


100.0 


6.14 


0.36 


100.0 


14.03 


0.53 


99.87 


8.36 


0.80 


99.83 


4.95 


0.45 



Measurement method abbreviations: (Standard) - Ground truth measurements based on three human observer measure- 
ments, (01-03) - Human Observers 1-3, (Gregson) - Gregson rectangle fitting UW, (HHFW) - Half Height Full Width [43] , 
(IDG) - ID Gaussian model fitting 14 4h (2DG) - 2D Gaussian model fitting \45\ , (ESP) - Extraction of Segment Profiles 
ffijl, (Graph) - Graph based method \24\ , (ARIA) - Autmated Retinal Image Analyzer [23] and (ETOS) - Edge Tracking on 
Orientation Scores. Dataset abbreviations: (KPIS) - the Kick Point Image Set, (CLRIS) - Central Light Reflex Imag e Set, 
(VDIS) - Vascular Disease Image Set and (HRIS) - the downsampled High Resolution Image Set (HRIS). See section \3.2.2\ 
for more details. 



and a scale too low only the small vessels). We found 
3 

that a = — 10 gave best results. For the CTOS algo- 
rithm we used a set of orientation scores constructed 
by Gabor wavelets at scale^jbased on r =5,10,15,20,25 

3r 

and 30, were the scale is given by a = — . 

Z7T 

Results of the tracking experiments are shown in 



Fig. 12 From this figure we see that, at complex sit- 
uations, the ETOS method (column 2 and 3) outper- 
forms the CTOS method (column 4). Best results are 
obtained when ETOS is used with invertible orienta- 
tion scores (column 2). The ETOS algorithm acting on 
the invertible orientation scores generated by the cake 
kernels only fails to correctly track blood vessel nr 5 
from image 3. The algorithm gives excellent results for 
all other vessels and manages to track the blood vessels 
through all complex situations, even when the contrast 
of the vessel edges is very low (as is at the case at the 
bifurcation involving seed point 3 in image 2). 

The performance of the ETOS algorithm is slightly 
decreased when applied to non- invertible orientation 
scores based on Gabor filters. It now fails to track 3 
vessels correctly: In test image 3, besides vessel nr 5, 
vessel nr 1 is incorrectly tracked and in test image 4 

4 The size of the Gabor wavelet is more intuitively de- 
scribed in terms of the spatial period r of the complex expo- 
nential. With ko = (k x , k y ) = (0, 3), the relation between the 



period r and the scale parameter a is given by a ■ 



k x T 



3r 
2tt' 



vessel nr 2 is not correctly tracked through the bifurca- 
tion. The scale selective property of the Gabor wavelets, 
resulting in non-invertible orientation scores, causes the 
ETOS algorithm to perform less accurate compared to 
the application to invertible orientation scores. 

The CTOS algorithm, which relies on a multi-scale 
orientation score approach, has the lowest performance. 
It fails to correctly track 5 blood vessels. In some cases, 
incorrect scale selection causes small parallel blood ves- 
sel to be detected as one large blood vessel. This effect 
is visible in test image 1 where blood vessel 2 and 4 and 
vessels 3 and 6 are merged at the end of their tracks. In 
test image 2 vessel segment nr 1 continues in vessel nr 
2 because of incorrect orientation detection. A similar 
problem occurs in image 4 where vessel nr 4 continues 
in vessel nr 2. 

In conclusion we can state that ETOS outperforms 
CTOS and that it gives best results when applied on 
invertible orientation scores. 

3.2.2 Validation of width measurements 

In the previous section we showed that the ETOS al- 
gorithm is very well capable of tracking blood vessels 
through complex situations. In this section we quanti- 
tatively validate the reliability of the measured vessel 
widths that are provided by the ETOS algorithm. This 
is done by comparing the measured widths to ground 
truth width measurements provided by the REVIEW 
database (6 . The REVIEW database consists of 16 
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color fundus images, which can be divided into 4 sub- 
sets: 1) Kick point image set (KPIS), 2) Central light 
reflex image set (CLRIS), 3) Vascular disease image set 
(VDIS) and 4) the high resolution image set (HRIS). 
Each image set represents images of different quality 
and resolution. Each image is provided with manual 
width measurements that are performed by three indi- 
vidual observers. A ground truth of the vessel widths is 
constructed by averaging the measurements of the ob- 
servers. The HRIS set contains high resolution images 
(3584x2438) and were downsampled by a factor of four 
before submission to the ETOS algorithm, such that the 
ground truth vessel widths are to be known with a sub- 
pixel accuracy of ±0.25 pixels. For more information on 
the dataset we would like to refer to [6 . When testing 
our algorithm, we tracked each segment by initializing 
the algorithm using the first pair of manually marked 
edge points. The same parameters that are described in 



section 3.2.1 were used. Fig. 14 shows a selection of the 
tracking results in comparison to ground truth vessel 
edge labeling. 

In total 5066 vessel width measurements are avail- 
able. The error between automated measurements and 
the ground truth measurements is defined as 



Xi = w i 



w? T = 



(35) 



where wi is the estimated width as measured by the 
ETOS algorithm (recall Eq. 21 ), and wf T is the ground 



truth width of the ith profile. To be able to compare 
our method with others we follow the same validation 
procedure as described in [6| l23[ l24] , where the main fo- 
cus is on the standard deviation of the errors. This is 
motivated by the idea that different implicit definitions 
of vessel widths may lead to consistent errors. If this 
bias however is consistent enough, the error could eas- 
ily be accounted for by subtraction of a bias constant. 
A low standard deviation of the errors indicates that 
the error is consistent. 

Table [T] shows the validation results of our ETOS al- 
gorithm, in comparison with methods by other authors 
that published their results using the same database [6l 
23 ,24]. The first four rows of the table show the results 
of the manual annotations (observer 1, 2 and 3) and the 
golden standard. The next four rows show results of four 
classic approaches to vessel width measurements: 

— Gregson: a rectangle is fitted to a vessel intensity 
profile, and the width is set so that the area under 
the rectangle is equal to the area under the profile 

E2]. 

- Half Height Full Width (HHFW): the standard half- 
height method, which uses thresholds set half-way 
between the maximum and minimum intensities at 
either side of an estimated center point [43] . 




10 15 

Ground truth width (pixels) 

Fig. 13: A scatter plot, plotting 5059 ground truth widths 
against the widths measured by our ETOS algorithm. The 
linear regression model y = 0.85 + 0.88x indicates an off- 
set of less then a pixel, suggesting that ETOS slightly over- 
estimates the vessel widths. The slope of 0.88 indicates a 
strong positive relation between the ground truth and mea- 
sured widths. 




Fig. 14: Several tracked vessel segments by ETOS in com- 
parison with manual width measurements. Left column shows 
the ground truth vessel edge labeling as provided by the RE- 
VIEW database. The middle column shows results obtained 
by the ETOS algorithm and the right column shows both the 
ground truth (in white) and our results (in red). 
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— ID Gaussian (IDG): a ID Gaussian model is fit to 
the vessel intensity profile [44] . 

— 2D Gaussian: a 2D Gaussian model is fit to the ves- 
sel intensity profile [45] . 

The next three rows give results of the most recent, 
state of the art methods that published their results: 

— The Extraction of Segment Profiles (ESP) is an ac- 
tive contour based algorithm developed by Al-Diri 
et al. 0. 

— The Graph method is a graph based edge segmen- 
tation technique developed by Xu et al. [24] . 

— The Automated Retinal Image Analyzer (ARIA) is 
an algorithm developed by Bankhead et al. [23] , 
they used a wavelet approach to vessel segmenta- 
tion after which the edge locations are refined. 

The last row shows the results we achieved using our 
ETOS algorithm. 

The column labeled with % shows the success rate, 
it indicates how many width measurements could suc- 
cessfully be validated. Each ground truth center point 
was matched with a center point generated by our vessel 
tracking algorithm. A pair of center points was regarded 
as successful when the distance between the two center 
points was smaller than half the ground truth vessel 
width and also if the detected point was not closer to 
another ground truth center point. The success percent- 
age is thus smaller then 100% whenever a measurement 
failed to converge. 

The column labeled with Mean indicates the mean 
vessel width of all the measured vessel profiles. The col- 
umn labeled with <j x indicates the standard deviation 
of the error (Eq. (K35|), a lower o~ x is favorable since it 



indicates that the error is consistent. 

From table [I] it can be observed that ESP, Graph, 
ARIA and our ETOS algorithm all outperform the clas- 
sic width measurement techniques. Also compared to 
the state of the art methods our algorithm scores very 
well. The ETOS algorithm performs remarkably well on 
the CLRIS dataset, which contains a large number of 
vessels with the central light reflex. For these images, 
the standard deviation of the errors is even lower than 
those of the observers. As for the other datasets, our 
method's performance is comparable to the state of the 
art methods. 



Fig. [13] shows a scatter plot, plotting the ground 
truth widths against the widths measured by our ETOS 
algorithm, together with a linear regression model that 
was fit through these points. The points are very much 
centered around the line y = x, indicating a strong 
positive correlation. This is confirmed by the slope of 
the linear regression model y = 0.85 + 0.88x, which is 
near to 1. The low number of outliers in the scatter 



plot confirms the low standard deviation in errors, as 
demonstrated by table [T] The offset of 0.85 indicates 
that the ETOS algorithm has the tendency to slightly 
overfestimate the vessel width. 

We conclude that our ETOS algorithm, which is 
highly capable of tracking blood vessels through all 
sorts of complex situations, also provides reliable width 
measurements. 



4 Vasculature tracking 

In this section we describe the extensions to the ETOS 
algorithm, which make it possible to construct mod- 
els of the complete retinal vasculature. Our vasculature 
tracking algorithm consists of: 

1. Optic disk detection. 

2. Seed point detection in the optic disk region. 

3. Correct Initialization of the ETOS algorith by ro- 
bust initial edge detection. 

4. Automatic termination based on a set of stopping 
criteria. 

5. Junction detection. 

6. Junction resolving. 

Each of these items is described in sections [4.1.1l to l4. 1.61 



and the complete algorithm is validated in section 4.2 



4.1 Methods 

4-1-1 Optic disk detection 

Since the blood vessels of interest all enter the eye 
through the optic disk, initial seed points can best be 
detected on a circle that is centered around the op- 
tic disk [rO][T2]. Our method therefore starts with op- 
tic disk detection. On color fundus images, the optic 
disk appears as a bright disk from which blood vessels 
emerge. The detection of the optic disk occurs in two 
phases. In the first phase a rough estimation of the optic 
disk position is made based on a method using variance 
filtering (Fig. [l5]o) ? as proposed by Sinthanayothin [46] . 
This method is based on the significant variance in pixel 
intensities that occur in the optic disk region. The esti- 
mated optic disk position is detected as the position of 
the global maximum in this image. For more detail on 
the implementation of this method we refer to [46] . 

Next, the estimated optic disk position is refined 
by edge focusing [47 on the optic disk boundary, fol- 
lowed by a weighted Hough transform [48] to find the 
parameters for a circle that best resembles the circum- 
ference of the optic disk. Here we assume that the optic 
disk is circular, and that it's size can be represented by 
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Fig. 15: Automated optic disk detection. Image (a) shows a partial view of a fundus image and (b) the corresponding variance 
image. The position of the global maximum of the variance image is used as a rough estimation of the optic disk center. In 
image (c) the red channel of a color fundus image around the estimated optic disk position is shown. The vessels in this 
sub-image are filtered out by a closing operator (d). To detect dominant edges, edge focussing is performed on the intensity 
profiles obtained form the lines depicted in figure (e). Dominant edge positions are shown as blue dots in image (f). A weighted 
Hough transform is performed using the these positions, the result is shown as a red circle. 



Edge detection, I v (r]) 




(a) 





(c) 



Fig. 16: Initial edge detection, (a) Based on local optima in the imaginary part of I 71 , potential vessel patches are formed 
(shown as horizontal blocks). The main patch of interest is denoted in red and neighbouring patches of interest in non- 
transparent blue, (b) The edges of the patches of interest are tracked in scale (lower part of this figure) up to the scale u where 
there are at least two potential edge pairs left (dashed line), (c) The strongest edges at this scale are initialized to be the vessel 
edges uq and vq. 



the radius of this circle. Edge focussing is performed 
on the red-channel for the reason that in this channel 
the blood vessels appear less dominant and the optic 
disk region appears brighter (Fig. [l5fc). To avoid dis- 
turbance of blood vessels during the detection of the 
optic disk boundary, the vessels are removed by a clos- 
ing operation with a circular structural element with 
a radius larger than the maximum vessel width, see 
Fig. p~5] l. Edge focussing is performed on a star-shaped 
set of profiles (Fig. [l5fc,f) and the detected edge posi- 
tions are used as input for the weighted Hough trans- 



form for circles. The weight of each edge position is 
determined by the scale up to which the edge can be 
traced in scale (edges disappear with increasing scale 
and the strongest edges survive longest when increas- 
ing scale in scale space). 

The optic disk radius Rod, found by the Hough 
transform for circles, describes the size of the optic disk. 
For adult human eyes, the average radius of the optic 
disk is known to be RHq^ = 0.92mm [49]. Vessels in the 
optic disk region typically have a caliber of (w)™f = 
0.15mm. While the resolution of retinal images may 
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vary from camera to camera, the physical dimensions of 
the human eye are rather constant between individuals. 
In order to make our algorithm applicable to all sorts of 
retinal images, we will parameterize our routines with 
a normalized distance measure 



(w)av = (w) 



Rod 1 



ref 

av re f 
n OD 



6 



Rod, 



(36) 



where, based on real physical values, (w) av describes 
typical retinal vessel calibers in pixels. We will use (w) av 
to make our vasculature tracking routines invariant to 
image resolution. 

4-1-2 Vessel likelihood map and seed point detection 

For the detection of the initial seed points, a vessel like- 
lihood map V : R 2 H> R of the optic disk region is 
constructed using invertible orientation scores. In an 
orientation score, isotropic structures in the image are 
uniformly distributed in the orientation direction of the 
score, thus providing a low response for all orientations. 
Elongated structures however are more localized in the 
orientation direction, meaning that an elongated struc- 
ture is mainly presented in the layers of the score which 
correspond to the structure's orientation. In the case of 
the color fundus images the elongated structures to be 
detected, the blood vessels, are dark structures com- 
pared to the background. Before images are subjected 
to our algorithms, they have their DC-component re- 
moved by means of high-pass filtering. As a result vessel 
pixels have negative values and background pixel values 
are assumed to be zero. In effect the real part of the ori- 
entation score at a certain vessel position and orienta- 
tion is negative, the orientation score is approximately 
zero at background area's. This is also demonstrated in 
Fig. |8b[ were the vessel in the score can be seen as a 
dark (negative value) blob on the plane. Based on these 
observations, we define a vessel likelihood map V(x) as: 



V(x) 



max Re( 

<96[0,2tt] 



-E//(x,0)). 



(37) 



Seed points are detected as local maxima on circu- 
lar intensity profiles with radii R = {RodA-5Rod}, 
where Rod is the detected optic disk radius, centered 



around the optic disk (Fig. 17a). A seed point is dis- 
carded whenever its value in the vessel likelihood map 
V is smaller then the average value of all points on 
the circle. For each remaining seed point Co, the ini- 
tial orientation Oo is detected as the orientation that 
provides the highest modulus of the orientation scores: 
#o = argmax|£//(c , 0)\. An additional filtering step, in 

6>G[0,2tt] 

which the seed points are classified as either true or 




(a) 



(b) 



Fig. 17: Initial seed point detection from vessel likelihood 
maps, (a) A vessel likelihood map of the optic disk region, 
with two circular profiles on which initial seed points are de- 
tected. Detected seed points are shown as red dots, discarded 
as black dots, (b) Edge initialization and true positive seed 
point selection. White arrows show the detected seed points, 
red arrows are seed points which are classified as false posi- 
tives. 



4-1-3 Initial 



detection 



false positive, is described in section 4.1.3 



The ETOS algorithm is initialized with a starting ves- 
sel center point Co, orientation Oo and edges uo and vo. 
Starting with an already detected initial center point Co 
and orientation O > intensity profile I 71 can be obtained 
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from the orientation scores using Eq. (25). Candidates 
edges are detected as local optima on the imaginary 
part of I 71 . Beside the main vessel edges that we are 
interested in, it is very likely that multiple other can- 
didate edges are detected as well (for example because 
of noise or a central light reflex). Therefore, we use an 
edge focussing approach to detect the dominant edges. 
Each combination of neighboring left and right edges 



will form potential vessel patches (see Fig. 16a). Note 
that a blood vessel with a central light reflex consists 
of two neighboring vessel patches. We start initial edge 
detection by detecting the main vessel patch of inter- 
est by scoring each edge pair (u,v), based on the initial 
center point estimate Co and initial orientation #o, as 
follows: 



l ||(u + v)/2-c | 
s CoA) (u,v) = ^ (u,v)e 2 (OMw)av) 2 



with 



1 f 1 

u,v) = - \U f (u + t(v-u),0)\dt. 

\\ u - v ll JO 



(38) 



(39) 



The function z^(u, v) provides the so called vessel value 
and basically is the average value of the modulus of 
the orientation score at orientation calculated from 
the left to the right edge. This value is high for elon- 
gated/vessel structures, and low for background struc- 



tures in the score. The exponential in Eq. (38) penalizes 



the distance of the edge pair to the initialized center 
point. The edge pair with the highest score will be our 
main pair of interest. 

Neighboring pairs are considered only if the distance 
between the nearest edges of the main and neighbor- 
ing pair is smaller than the width of the main patch, 



Fig. 16a shows typical results. All edges from the main 
patch and its neighbors are now traced in scale by gaus- 
sian blurring the profile I 71 . The edges are traced up to 
the scale where there are still at least two possible edge 



pairs left (dashed line Fig. 16b). The strongest edges at 
this scale are chosen as the true vessel edge points Uo 
and vo (Fig. |16c[ ). 

An initialized seed point whose vessel value is lower 
than a threshold T v ^ is regarded as a false positive 
(Fig. 17b). The threshold T v ^ is defined as: 



(40) 



where using all Nb initialized branches, each with ini- 
tial edge points and orientation (uq, Vq and 0q resp.), 
an average vessel value (v) av is calculated, and j3 v is a 
sensitivity parameter (typically between and 1.5). All 



other seed points are regarded as true positive and are 
submitted to the ETOS algorithm to start expanding 
a model of the retinal vasculature. Note that a higher 
value of j3 v allows more seed points to be submitted to 
the ETOS algorithm. 

4-1-4 Stopping criteria 

For ETOS to stop tracking a single vessel, three stopping- 
criteria are defined: 

1. Tracking stops whenever a vessel being tracked leaves 
a prescribed region of interest. For this purpose a 
mask image is generated that covers the (typically 
circular) field of view in the fundus image. The mask 
is eroded by 40 pixels to eliminate potential border 
effects. During each iteration k it is checked whether 
or not the center point lies within the mask; if 
not, the algorithm stops. 

2. Tracking stops whenever a blood vessel is already 
tracked. Each tracked segment is used to construct 
a pixel map, in which each pixel within the vessel 
edge contours is set to 1 and outside to 0. Whenever 
a point within the vessel being tracked lies within 
the pixel map for \4(w) av /X] iterations in a row, the 
ETOS algorithm terminates. Recall (w) av defined in 
Eq. (36 and A being the step size (Fig. [6]). 

3. Tracking stops whenever the vessel value vq (Eq. ([39]) ) 
drops below a certain threshold value T Uja : 



T U:a = a v {v) c 



(41) 



where ct v is a sensitivity parameter (typically be- 
tween and 1.5). Here the assumption is made that 
the average vessel value (v) av , calculated in the op- 
tic disk region, is a good indicator of the vessel val- 
ues of the vessels elsewhere in the retina. Increasing 
a v may cause the algorithm to stop too early, low- 
ering the value will increase the risk of continuing 
vessel tracking into the background). 

4-1-5 Junction detection 

In order to model the complete vasculature, starting 
by expanding a model from a set of initial seed points, 
the vasculature tracking algorithm should also be able 
to automatically detect junctions. Junction points are 
points where blood vessels bifurcate/branch or where 
two blood vessels cross. Either way, at a junction point 
multiple orientations may be observed. At a point on a 
straight vessel the modulus of the orientation column 
ideally contains two local maxima; one in the positive 
direction, + , and one (with a tt phase difference) in 
the negative direction, 6_. A candidate vessel junction 
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Fig. 18: Junction detection. The top image: A tracked ves- 
sel segmented in white with candidate junction points indi- 
cated by colored arrows. The junction points are clustered on 
position (indicated by white ellipsoids) after which they are 
clustered on orientation (indicated by different colors). Bot- 
tom image: Each cluster is merged to a single junction point 
and based on proximity to other junction points, orientation 
and width, junction points are classified as bifurcation ( green 
arrows) or crossing (red dashed arrows). 



point at the left edge is detected whenever there exists 
another local maximum in between 9+ and 0_. Simi- 
larly, a candidate junction point is detected at the right 
vessel edge whenever another local maximum is present 
between #_ and Q+. 

During vessel tracking, the orientation columns at 
the left and right edge are scanned for the presence of a 
junction point. If a junction point is found, the location 
and orientation are stored, see the top image of Fig. 18 
The detected candidate junction points are clustered on 
position by grouping all points whose distance to one 
another is smaller than (w) av . Within each cluster, the 
candidate junction points are clustered on orientation 
to prevent merging two proximate junction points (e.g. 
a branching point near a crossing point). Clustering on 
orientation is done according to the number of local 
maxima in the histogram of orientations. Finally each 
separate cluster is merged to a single junction point by 
averaging the positions and by taking the most com- 



mon orientation within the cluster. The found center 
points and orientations are then subjected to the edge 
initialization method described in section [4.1.3[ and are 
discarded whenever their vessel value is lower than T v $. 

Each junction point is classified as a crossing point 
or bifurcation point. A crossing point is classified as 
such when there are two junction points on each side 
of the blood vessel, which are closer to each other then 
two times the width of the vessel, and of which the 

7T 

orientation difference is in the range tt ± — . A typical 
result is shown in the bottom image of in Fig. [18| 

Each detected seed point is assigned two ID num- 
bers, a number that is unique for each vessel segment 
and the ID number of the vessel from which it origi- 
nates. This way the relation between each vessel seg- 
ment remains known, and the vessel segments can be 
organized in a hierarchical fashion. Each detected bi- 
furcation point is stored in a list of untracked seed 
points, each detected crossing point is stored in a list 
of untracked crossing-points. Whenever a seed point is 
tracked, it is removed from the list. When all the bi- 
furcations are evaluated, the crossing points are con- 
sidered. A crossing point is discarded if it is already 
tracked, it is tracked otherwise. Vasculature tracking 
terminates whenever all detected bifurcations and cross- 
ings are evaluated. 

4.1.6 Junction resolving 



As mentioned in section |4.1.4| vessel tracking is termi- 
nated whenever the algorithm is tracking a vessel that 
is already tracked. This criterion can be met at several 
situations, where for each situation appropriate actions 
need to be taken in order to maintain correct topolog- 
ical models of the vasculature. A detected point that 
suggests inappropriate modeling of the vasculature will 
be called an unresolved junction point. The detection of 
an unresolved junction point, together with the corre- 
sponding actions that are necessary to solve it will be 
called junction resolving. 

The junction resolving steps are described using the 
following nomenclature: S Q id is the tracked segment of 
the already established track before the junction point 
(Source), T u is the tracked segment of the same track 
after the junction point (Target), S new is the segment 
of the new track before the junction point, and T new is 
the segment of the new track after the junction point. 
Our junction resolving algorithm is able to recognize 
nine different types of unresolved junction points based 
on a sequence of three questions: 

1. Where on the original (old) track and where on the 
new track does the unresolved junction point occur? 
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(a) 



(b) 



(c) 



Fig. 19: An example of the detection of, and solution to, unresolved junction points, (a) Two tracks end up tracking the same 
blood vessel, resulting in an unresolved junction being identified, (b) Based on vessel width, orientation and pixel intensities 
the correct track (red dashed) is identified. As an immediate action the incorrect track will be cut-off at the junction point and 
is marked as being open-ended. During the postponed actions, when all vessels are being tracked, the open-ended track (white) 
is matched with a vessel that starts in the extend of it (blue). The result of the junction resolving algorithm is shown in (c). 



2. Are the two target tracks different to each other, and 
if so, which one is correct? For example, two tracks 
can overlap while being different in vessel width if 
one of the tracks was incorrectly initialized. 

3. If the two target tracks are not different to each 
other, which one of the two source tracks belongs to 
the target track? 

More detail on the implementation and examples of 
the junction resolving algorithm can be found in [50] . 
Within this paper we only briefly demonstrate the al- 
gorithm with one example. 



In our example, illustrated in Fig. 19 , the ETOS al- 
gorithm fails to correctly track a blood vessel through a 



crossing (white track in Fig. |19a). Now when the other 



vessel is being tracked (red dashed track in Fig. 19a), 
our algorithm will detect that from a certain point two 
tracks cover the same vessel. This point will be recog- 
nized as an unresolved junction point. To answer the 
questions, starting with the first: The unresolved junc- 
tion point occurs somewhere halfway on the already 
established track, meaning we have source and target 
vessel segments (S i d , T oM , S new and T new ). Answering 
question two: The segments after the junction points, 
the targets T u and T newi represent the same blood 
vessel. Since the tracks before the junction point (the 
sources) do not overlap, one of the tracks is incorrect. 
Answering question three: Based on vessel width, orien- 
tation and pixel intensities, the target track is matched 
to the source of the new track. 

Now that the problem is identified, appropriate ac- 
tions need to be taken. The actions to solve unresolved 
junction points can be categorized in immediate ac- 
tions and postponed actions. The immediate actions are 



taken immediately after the unresolved junction point 
is identified and they involve, among others, deleting 
incorrect /abundant (parts of) vessel models and merg- 
ing vessel models. The postponed actions are postponed 
until all blood vessels are tracked and they involve, 
among others, merging vessel models and reassigning 
the parent IDs of the vessels of which the parent ves- 
sel was deleted during one of the immediate actions. 
A complete overview of the immediate and postponed 
actions can be found in [50 j. In our example, the im- 
mediate actions involve merging the source S new of the 
new track with the part of the already established track 
after the junction point, T ^. The old track now only 
consists of S id (see Fig. 19b) and is marked as open- 



ended because it is abruptly cut of at a potential cross- 
ing point. During the postponed actions, tracks marked 
as being open-ended are matched with tracks that start 
in the extend of the cut-off tracks, in Fig. |19b| such a 
potential match is denoted with T. If a good match can 
be found, based on proximity, vessel width and orien- 



tation, the two tracks are merged (Fig. 19c) 



4.2 Validation 



In section 3.2.2 we demonstrated the reliability of the 
width measurements provided by the ETOS algorithm. 
In the following section we validate the topological cor- 
rectness of the complete vasculature models that are 
generated by our algorithm. The correctness of the mod- 
els is validated by analyzing the junction points. The 
results discussed in this section are generated with the 
same parameters for the ETOS algorithm that are de- 



scribed in section 3.2.1 The vessel value sensitivity pa- 
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(a) Vessel tree extraction (b) Junction points (c) Junction points 



Fig. 20: (a):The hierarchical structure of the generated vasculature models allow the segmentation and analysis of complete 
branches, (b) The automatic extraction of branching (yellow) and crossing points (red), (c) The distance to the optic disk; 
a feature that can easily be extracted because of the guaranteed connectedness of vessel segments in the generated vasculature 
models. 



Table 2: Validation of detected bifurcations 



Disease 




Correct 


Bifurcations 

Part of a crossing 


False vessel 


Total 




Crossings 

Correct Incorrect 


Total 


Healthy 


90 


(81.82%) 


16 (14.55%) 


4 (3.64%) 


110 


41 


(100.00%) 


(0.00%) 


41 


Diabetes 


111 


(73.03%) 


25 (16.45%) 


16 (10.53%) 


152 


37 


(88.10%) 


5 (11.90%) 


42 


Glaucoma 


89 


(74.79%) 


30 (25.21%) 


(0.00%) 


119 


31 


(100.00%) 


(0.00%) 


31 


All 


290 


(76.12%) 


71 (18.64%) 


20 (5.25%) 


381 


109 


(95.61%) 


5 (4.39%) 


114 



rameters a vi Eq. (41), and /3 U , Eq. (40), are both set to 



0.5. A typical model generated by the vasculature track- 
ing algorithm with these parameters, together with po- 



tential applications is shown in figure 20 



4.2.1 Validation of topological correctness 

The hierarchical structure of the created vasculature 
model allows easy extraction of junction points. For 
each vessel it is known from which parent vessel it orig- 
inates; bifurcations can thus be directly extracted from 
the model. Crossings can easily be extracted by detect- 
ing overlapping vessel segments. Fig. |20b| provides an 
overview of detected junction points for one image of 



the HRFI-database, Fig. [21] shows detailed views of a 
bifurcations and a crossing. We evaluated the junction 
points of the first three images from each of the three 
HRFI-datasets (healthy, diabetes and glaucoma). The 
9 generated vasculature models provided 495 junction 
points, of which 381 were bifurcation points and 114 
were crossing points. Bifurcations were marked as in- 
correct whenever they were actually part of a crossing, 
or whenever the vessel originating from the bifurcation 
did not represent a blood vessel according to the ground 
truth pixel map provided by the HRFI-database. This 
last type of error indicates the presence of inappropri- 
ate single vessel models, i.e. models that do not repre- 
sent true vessels but rather noise, non-vessel elongated 
structures such as the optic disk border or pathological 



features such as aneurysms. Crossings are extracted by 
searching the vasculature model for overlapping vessel 
segments. A false positive crossing can thus only occur 
if false positive vessel segments exists within the mod- 
els. The results are summarized in table |2j 



In total 290 bifurcations and 109 crossings were cor- 
rectly detected, corresponding to precision rates of 76.54% 
and 96.03% respectively. Most of the incorrectly iden- 
tified bifurcations were correct in the sense that they 
represent a true blood vessel, however they were actu- 
ally part of a crossing. Only 5.25% of the bifurcations 
were incorrect in the sense that they did not represent 
a blood vessel. While the first kind of false positive bi- 
furcations only affect the topological correctness of the 
model, the latter also pollutes the model with false pos- 
itive vessel segments. The low percentage of false pos- 
itive bifurcations that did not represent a blood ves- 
sel indicates that the generated vasculature models are 
very clean in the sense that almost all vessel segments 
actually represent true blood vessels. The false posi- 
tive bifurcations which were actually part of a cross- 
ing, generally arise whenever the source of the crossing 
vessel is not detected. This might for example occur if 
the vessel enters the fundus image at the image bound- 
ary. As described in section [4. 1.5[ even if two junction 
points are correctly identified as a crossing, it is later re- 
identified as two bifurcations if the crossing vessel still 
is not tracked after all bifurcation points are tracked. 
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This then leads to the inclusion of two false-positive 
bifurcation points in the model. 




Fig. 21: A typical bifurcation (left) and crossing (right), de- 
tected from a model generated by our vasculature tracking 
algorithm. 



5 Conclusion and future work 

Inspired by the perceptual organization of orientations 
in the early visual system of higher mammals, we pro- 
posed two novel vessel tracking methods based on ori- 
entation scores: ETOS and CTOS. In this paper we 
demonstrated that by representing image information 
based on position and orientation in an orientation score, 
one can exploit the disentanglement of crossing struc- 
tures to track blood vessels through crossing points. 
The ETOS algorithm tracks vessel edges through the 
orientation score of an image and can generally be used 
with both invertible and non-invert ible orientation scores, 
which were in this paper constructed with cake wavelets 
and Gabor wavelets respectively. We demonstrated that 
best results were obtained using invertible orientation 
scores. The CTOS algorithm is a centerline tracking 
method and is especially designed to work on a multi- 
scale set of Gabor wavelet based orientation scores. 
While this algorithm is very fast, the multi-scale ap- 
proach makes the algorithm less stable at critical ves- 
sel points (crossings and parallel vessels) compared to 
ETOS. 

The ETOS algorithm was used as a basis for our 
vasculature tracking algorithm, with which we can con- 
struct detailed hierarchical models of the retinal vas- 
culature. Information about the vessel orientation and 
width, as well as the location and geometry of bifurca- 
tion and crossing points can directly be extracted from 
the models. Within this paper we validated the reliabil- 
ity of the width measurements provided by the models 
using ground truth data, and showed that our method 
performs excellent in comparison with other state of 
the art algorithms. Validation of the topology of the 



models showed that our method constructs clean mod- 
els, ie. they contain very few false positive vessels. The 
validation also showed that the main source of incor- 
rect topological relations was the miss-classification of 
crossings as two separate bifurcations. As part of fu- 
ture research we will therefore focus on increasing the 
topological correctness of the generated models both by 
introducing "post-tracking" analysis methods, analysis 
and correction of the models after vasculature track- 
ing is done, and improving the existing " intra-tracking" 
analysis methods such as stopping criteria and junc- 
tion resolving algorithms. On the theoretical side we 
will continue exploring active contour models within 
(enhanced [40]) orientation scores with application to 
retinal images. 

A A mathematical underpinning of 
optimization in the e^ — e# tangent planes V 

Before we provide the key-geometrical principle behind our 
ETOS algorithm, we make some comments on the moving 
frame of reference that lives in the domain of an orientation 
score. 

The domain of an orientation requires a curved (non- 
Euclidean) geometry as can be seen in Figure ^ To this end 
we identify the coupled space of positions and orientations 
M 2 xi S 1 with SE(2). This needs a bit of explanation. To see 
whether a local line element (x',0') G M 2 x S 1 with posi- 
tion x / E M 2 and orientation f is aligned with its neigh- 
bors one applies a rigid body motion g = (x,R#) G SE(2), 
with counter-clockwise rotation and translation via Eq. (m). 
This puts a 1-to-l correspondence between rigid body mo- 
tions (x,y,Ho) and elements from I 2 x 5 1 as every local 
line element (x, y, 6) can be obtained by a local line element 
(0,0,0) positioned at the origin = (0,0) pointing in say 
x-direction: 

(x, 6) = (x, R e )(0, 0) ^ (x, 6) = (x, R e ). 

If we do so we see that the curved domain of an orientation 
score is due to the non-commutative group structure of the 
rigid body motion (also known as roto-translation group or 
Euclidean motion group) SE{2) in the plane. Indeed, a con- 
catenation of two rigid body motions is again a rigid body 
motion and this induces the group-product given by Eq. (|8| 
It is a semi- direct product of rotations SO (2) and trans- 
lations M 2 , as the rotation part affects the position part. As a 
result this group product is non-commutative, i.e. gg' can be 
very different from g' g. Therefore it is conceptually wrong to 
consider the domain of an orientation scores as the flat Carte- 
sian space I 2 x 5 1 with Euclidean metric and it is wrong to 
just take derivatives of orientation scores only in d x direction 
or only in d y direction. For formal underpinning of this state- 
ment see [2lJ Thm.21], for a short illustration see [40j Fig. 2. 6]. 
Instead one must differentiate along the moving frame of ref- 
erence {de , d% , d v } with 

de,d^ = cos 0d x + sin 0d y , d v = — sin 6d x + cos 0d y . 

This a priori moving frame of reference coincides with the so- 
called Lie algebra of left-invariant vector fields on the group 
SE{2). The curved geometry in the orientation score requires 
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us to apply parallel transport along this moving frame of 
reference. For instance in Figure 1(d) we see examples of auto- 
parallel curves whose tangent vector is always constant to the 
moving frame of reference. The non-commutative structure of 
SE(2) (due to the fact that rotations and translations do not 
commute) is reflected by the commutators in the Lie-algebra 

[d e , d£] = d v , [de,d v ] = -dz , [d£,dr,] = 0, 

where [A, B] = AB — BA. As a result first moving in £ and 
then moving in 6 direction brings us to a different point in 
SE(2) = R 2 x S 1 than first moving in direction and then 
moving in £ direction (as can be deduced in Figure 1(b)). 



A.l The geometrical principle behind the 
ETOS-algorithm 



The ETOS-algorithm presented in section |3.1.1| heavily re- 
lies on local optimization in each 2D-tangent plan^J V = 
span{<9#, c^} spanned by de and d v = — sin0d x + cos0d y in 
tangent-bundle 

T(SE(2)) = |J T g (SE(2)), 

geSE(2) 

where T g (SE(2)) denotes the tangent space at g = (x, y, 6) G 
SE{2). To this end we recall Figure [5] where one can observe 
that each tangent vector to the lifted curve s i— >> (x(s), y(s), 0(s)) 
with 6(s) — arg(£ (s)-\-iy(s)) (e.g. the curve following the right 
edge of a blood vessel) in SE{2) is pointing orthogonal to a 
plane V (plotted in yellow). 

Within V we see that the maximum values of the absolute 
value of the imaginary part of the score is located at the origin 
of the yellow plane in the tangent space. 

Definition 1 A smooth curve s \-> j(s) — (x(s),y(s),6(s)) 
in SE(2) is called the lifted curve of a smooth planar curve 
iff for all s G [0,£] we have 6(s) = arg(x(s) + iy(s)). 

Definition 2 If a curve 7 is equal to the lift of its spatial 
projection (i.e. if its satisfies Eq. (|42|)) it is called horizontal. 



Tangent vectors to horizontal curves always lay in the tangent 
plane H = span{<9#, := cos0d x + sin0d y }, spanned by d% 
and da, since one has 



6{s) = arg(x(s) + iy(s)) o 
7(s) e span{a c | 7(s) ,d e }, 



(42) 



In this appendix we will underpin and discuss our funda- 
mental venture point: The most salient curves in the smooth 
imaginary part/real part /absolute value C : SE{2) — > R of 
the smooth orientation score U : SE{2) — >• C are given by 

7 is horizontal such that deC(j) = and (d v \ C)(pf) — 

The intuitive idea behind this is as follows. Let C : SE{2) 1— >• 
R + denote an a priori given cost. Now, if a horizontal curve 
s 1 — V 7(s) := (x(s),^(s)) satisfies 

d v C(yi(s),9(s)) = (- sinO(s)d x C + cos6d y C)(x(s),6(s)) 
= d C(pc(s),G(s)) = O 

5 tangent vectors can be considered as differential operators 
acting on smooth locally defined functions [51] . In our case 
this boils down to replacing e^ by d x , e y by d y and e# by 



(43) 

with ^C(x(s), 6{s)) < and 0gC(x(s), 0(s)) < 0, then there 
is no gain in moving] the curve s \-> 7(5) = (x(s),0(s)) 
in directions orthogonal to the spatial propagation vector 



ft 



'ehoo 



: cos 0{s)d x + sin 6{s)d y 



In fact, we expect the curve 7(5) = (x(s),0(s)) satisfying 
( |43] ) to be optimal in some sense within the sub-Riemannian 
manifold 

M = (S'£;(2),span{a c ,%},C7 / 3) 

with metric tensor Gp : SE(2) x T(SE{2)) x (T(SE{2))) R 
is given by 



(x, y ,e) 



dO d(9+ 

f3 2 (cos 6 dx-\- sin 6 dy) ® (cos 0dx-\-s'm Qdy), 



(44) 



where denotes the tensor product. Note that 



(cos #dx + sin 0dy)(7(s)) = x(s) cos 0(s) + y(s) sin 6(s) ^ 
Gp(j(s),j(s)) = f3 2 \ cos 0(s)x(s)+sine(s)y(s)\ 2 + \9(s)\ 2 . 

If s equals arclength (i.e. we have ||x(s)|| 2 = 1) of the spatial 
part s 1— >- x(s) = (x(s),y(s)) of the horizontal curve s \-t 
j(s) = (x(s),0(s)) we have k 2 (s) = |^(s)| 2 , so that 



ll-VWII = 7 G /3(7(*),7(s)) = y/K?(8) + l3*. 



(45) 



Intuitively, such a sub-Riemannian manifold M equals the 
group SE{2) where one restricts oneself to horizontal curves 
with a constant relative penalty for bending and stretching 
determined by f3 > which has physical dimension one over 
length. 

For special cases of C we can show that our geometrical 
principle indeed produces optimal curves in SE(2). 



A. 2 Application of the geometrical principle to 
completion fields 

In general the real part, imaginary part, or absolute value 
of an orientation score is a complicated function on SE{2). 
Therefore we will consider the case where C : SE(2) — > R+ 
is a so-called "completion field" 52 , 53 , 54 , 55 , 56 . This cor- 
responds to collision probability densities of a single source 
particle gi E SE(2) and a single sink g 2 G SE(2). 

There exist remarkable relations 39 , 57 between optimal 
curves (i.e. curves minimizing an optimal control problem on 
SE(2)) and solutions of Eq. ( |43| ) for special cases where C 
denotes a so-called completion distribution (or "completion 
field"). Given two sources at the origin (0,0,0) and at gi — 
(xi, yi,0i), such completion fields are defined as products of 
resolvent Green's functions of stochastic processes for contour 
completion |58] and contour enhancement 59,60 in SE(2): 



C{g) := A 2 • Rx(9r 1 9)RxA92 1 9), 



(46) 



where R\(g) denotes the probability density of finding a ran- 
dom walker g in the underlying stochastic process [59,58 
given that it started at g = (x,y,0) — (0,0,0) regardless 
its memoryless traveling time T which is negatively exponen- 
tially distributed with expectation E(T) = A -1 . Furthermore, 

6 In such a way that the perturbed curve is again horizontal. 
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Contour completion 





Contour enhancement 





Fig. 22: Top: A contour completion process and the cor- 
responding Green's function (in SE{2)) and its planar 0- 
integrated version. Bottom: A contour enhancement process 
the corresponding Green's function (in SE{2)) and its pla- 
nar 6-integrated version. In the completion process one has 
randomness in with variance 2Dn > and non-random 
advection in £. In the enhancement process one has random- 
ness both in 6 -direction (with variance 2Dn > 0) and in 
^-direction (with variance 2L>22 > 0). 



g M> R\^(g) denotes the adjoint resolvent kernel (i.e. the re- 
solvent that arises by taking the adjoint of the generator [571 
ch:4.4]. 

For a concise overview on (Green's functions of ) contour 
completion and enhancement see |61j . Exact formulas for the 
resolvent Green's function R\ for contour enhancement can 
be found in [59], whereas exact formulae for resolvent kernels 
R\ for contour completion (direction process) can be found 
in [57] . In both cases there are representations involving 4 
Mathieu functions. For a visual impression of exact Green's 
functions see Figure [22] (where in dashed lines we have de- 
picted level sets of the corresponding Heisenberg approxima- 
tions that we will discuss and employ in the next subsection). 

These completion fields relate to the well-known Brown- 
ian bridges (in probability theory) where the traveling time 
is integrated out. This relation is relevant, since it is known 
that such Brownian bridges concentrate on geodesies |62| . 

If G t : SE(2) — » R+ denotes the time dependent Green's 
function of the Fokker-Planck equation of the underlying time 
dependent stochastic process [6T1I591I58] at time t > we have 

oo t 

C(g) = X 2 J J G t - s (g- 1 g)e-^ t -^G s (g- 1 g)e- Xs dsdt . 
o o 

This identity follows from two facts. Firstly, the resolvent 
Green's function follow from the time dependent Green's func- 
tion via Laplace transform with respect to time 

R x (g) = \C(t^G t (g))(\). 



Secondly, a temporal convolution relates to a product in the 
Laplace domain. As Brownian bridge measures concentrate 
on geodesies when A — » 0, cf. [62], [63] App.B], the com- 
pletion field for contour enhancement concentrates on sub- 
Riemannian geodesies (shortest distance curves) within M. 

Such a sub-Riemannian geodesic 7 = (x, 9) connecting g\ 
and #2 (with g\ and g<z sufficiently close^J) is the lifted curve 
of the minimizer to the following optimal control problem 



min f Jk 2 (s)+(3 2 ds 

xGCHlO,^! 2 ), 

£ > 0, 
(x(0),x(0))=flfi, 
(x(£),x(£))=02 



(47) 



7 =(x,^)eC 1 ([0,45E(2)), 
£ > 0, 
(x(0),x(0))=si, 

(XW,XW)^2, 

0(s) = arg(x(s) + iy(s)) 

with spatial arc-length parameter s > and with free total 
length £ > and with curvature k,(s) = \\x(s)\\ and where 
metric tensor G^ is defined by Eq. (|44j). 

On the other hand, in his paper 58 Mumford showed that 
the modes of the direction process (also known as the contour 
completion process) coincide with elastica curves which are 
the solutions to the following optimal control problem 



inf 

x g C^O, £],£ > 
(x(0),x(0)) = <? 0; 
(x(L),x(*))=pi 



1 



(s) + /3 2 ds 



(48) 



Similar to the Onsager-Machlup approach to optimal paths 
|67j he obtains these modes by looking at the most proba- 
ble/likely realization of discretized versions of the direction 
process. As this cannot be (efficiently) realized in practice, 
one needs a more tangible description of the mode. To this 
end, we will call solution curves of ( |43| ) the "modes" . In case of 
the direction process and in case one uses the completion dis- 
tribution (Eq. ( |46[ )) for the function C(g), the solution curves 
of ( |43| ) indeed seem to numerically coincide with the elastica. 
This experimentally underpins our conjecture in [57], where 
we have shown such a result does hold for the corresponding 
Heisenberg approximations as we will explain next. 



A. 3 In the Heisenberg approximation of completion 
fields our approach produces B-splines 

The Heisenberg group approximation (obtained by contrac- 
tion 59 ) of the Green's functions and induced completion 
field arises by replacing the moving frame of left-invariant 
vector fields 

{cos Odx + sin 6d y , — sin 6d x + cos 0d y , do } 

on SE{2) by the moving frame of reference of left-invariant 
vector fields on the Heisenberg group 



{d x + 0d y , d y ,de} 



7 This means g^ 1 gi is within the range of the exponential 
map of the sub-Riemannian control problem, i.e. g\ and g^ 
can be connected via a stationary curve without cusps |641 



mm. 
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and by replacing spatial arc-lengt parametrization via s by 
spatial coordinate x. Intuitively, such replacement boils down 
to replacing the space of positions and orientations by the 
space of positions and velo citi es. 

When contracting Eq. ( |46| ) towards the Heisenberg group 
H3 we obtain 



ay dO 



(49) 



where again in the Heisenberg group each tangent plane {do , d y } 
is orthogonal to the propagation direction d x -\-0d y and where 
the completion field is the product of two resolvent Green's 
functions 



C(x,y,6) = R x (x,y,6)R x (-x + x 1 ,y- 
which can be derived in exact form as 
AV3 



i(x - xi), -6) 



Rx(x,y,6) 



2D 1i7 tx 2 



_ 3(x0-2y)^+x 
4*3/?!! 



u(x), 



where Dn > stands for the amount of diffusion in 6- 
direction and u for the unit step function, for det ails^] see 
(cf. 57, Thm 4.6]). Interestingly, the solution of (|49| is a 
third order polynomial y(x) naturally lifted to (the corre- 
sponding sub-Riemannian manifold within) H{3) by setting 
9(x) — y'(x). The solutions are therefore cubic B-splines 
which are the solutions of the Euler-Lagrange equation 

y^\x) = 6^\x) 







of a curve optimization problem which arises by contracting 
( |48| ) towards the Heisenberg group. This gives the following 
"Heisenberg group equivalent" of (148}: 



7 (-) = (-,2/(-)^(-))€C 1 (H 3 ) 
7(0) = (0,0,0), 
7(012) = (22, 2/2, #2), 
6{x) = y'(x), 



X 2 

IP 2 





+ \Q'(x)\ 2 dx 



4(3y 2 2 +3x 2 y 2 2 +x 2 2 6 2 2 ) 



which takes the minimum along a lifted cubic-B spline 
(x, y(x), y' (x)) (with y{x) a third order polynomial matching 
the boundary conditions). For details see 57, Eq. 4j62] and 
[631 ch:9.1.1]. See Figure [23] 



A. 4 Concluding remarks 

By the results of the previous section the con jecture rises 
whether elastica curves coincide with Eq. \4Q\ as such re- 
lation holds for their counterparts in the Heisenberg group. 
Numeric computations seem to provide a confirmation of this 
conjecture, see Figure [24] 

On the one hand, we expect with respect to the contour 
enhancement process that our approach produces the sub- 
Riemannian geodesies (based on the results in [62ll63l[39] ), 
but this a point for future investigation. On the other hand, 
the conjecture together with the result in [58] (and the result 
that B-splines solve Eq. \4Q\ would underpin our altern ativ e 
applicable definition of modes as solution curves of Eq. ( |43| ). 

In fact this means that optimization in each (77, #)-plane 
V , as is done in our ETOS algorithm, produces the most prob- 
able curves in direction processes (for contour completion) in 
the sense of [581167] ! 




Fig. 23: The intersection of the planes {(x,y,6) G 
R 3 I d e C(x,y,0) = 0} and {(x,y,Q) G M 3 | d y C(x,y,0) = 0} 
of the Heisenberg a ppro ximation C(g) of the completion field 
C(g) given by Eq. Ujfy, produces a cubic B-spline lifted in 
H(3), i.e. (x, y(x), 0(x) = y' (x)) withy^(x) = 0. Boundary 
conditions have been set to x\ — yi = 0, 61 = OA, X2 — 2, 
yi = 0, 9 2 = -OA, D 1± = 1/8. 




d v C — 


1 \ 

9 




V 




8 Set = = in [57]. 



Fig. 24: The intersection of the planes {g G 
SE{2) I deC(g) = 0} (depicted in red) and 
{g G SE{2) I d v C(g) = 0} (depicted in yellow) of the 
exac t co mpletion field g = (x,y,6) \-> C(g) given by 
Eq. \4ty ) with xi = y± = 61 = and X2 = 3, y2 = —1, 
9 2 = 7/4vr A = 0.1, Dn = 1/32. The intersection of 
these planes seems to coincides with an elastica (with 
f3 2 — ^\D\\), which we plotted in dashed red in the top 
figure and by white balls in the bottom figures 
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